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We study synchronization of sinusoidally coupled pliase oscillators on networks with modular structure and 
a large number of oscillators in each community. Of particular interest is the hierarchy of local and global 
synchrony, i.e., synchrony within and between communities, respectively. Using the recent ansatz of Ott and 
Antonsen, we find that the degree of local synchrony can be determined from a set of coupled low-dimensional 
equations. If the number of communities in the network is large, a low-dimensional description of global 
synchrony can be also found. Using these results, we study bifurcations between different types of synchrony. 
We find that, depending on the relative strength of local and global coupling, the transition to synchrony in the 
network can be mediated by local or global effects. 
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I. INTRODUCTION 

Large networks of coupled oscillators are pervasive in sci- 
ence and nature and serve as an important model for studying 
emergent collective behavior. Some examples include syn- 
chronized flashing of fireflies yj, cardiac pacemaker cells 
1 2], walker-induced oscillations of some pedestrian bridges 
1 3], Josephson junction circuits 101, and circadian rhythms 
in mammals isj]. A paradigmatic model of the emergence of 
synchrony in systems of coupled oscillators is the Kuramoto 
model [|6[|, in which each oscillator is described by a phase 
angle 0„ that evolves as 
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where a;„ is the intrinsic frequency of oscillator n, Anm rep- 
resents the strength of the coupling from oscillator m to os- 
cillator n, and n,m = 1, . . . , A^. The classical all-to-all Ku- 
ramoto model corresponds to Anm = k. The study of gen- 
eralizations of the Kuramoto model has become an important 
area of research. Some examples of such generalizations in- 
clude systems with time-delays |01, network structure 
non-local coupling ifioll . external forcing ifTlll . non-sinusoidal 
coupling lfl2ll . cluster synchrony ifTsll . coupled excitable os- 
cillators uM, bimodal distributions of oscillator frequencies 
ifTsll . phase resetting ||T6ll . time-dependent connect ! vitv f lTll. 
noise jTstl , and communities of coupled oscillators lll9l - l22ll . 

In this paper we study the case where the coupling strength 
is not uniform, but rather defines a network that has strong 
modular, or community, structure. Synchrony on heteroge- 
neous networks has been studied in the past, both for phase os- 
cillator systems [8] and other dynamical systems [231. Much 
recent work has focused on the synchronization of phase os- 
cillators on networks with modular structure |[l9l - l22ll . While 
the link between community topology and synchronization is 
well established ll24ll . there are few analytical results that de- 
scribe synchronization in modular networks. Reference ifloll 
developed a framework to study a general number of commu- 
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nities, assuming that oscillators within communities are iden- 
tical. Reference [20] analyzed the linear stability of the inco- 
herent state for a system of coupled communities of hetero- 
geneous phase oscillators. The same system was considered 
in Ref. [25], where a set of coupled low-dimensional equa- 
tions governing the dynamics of the community order param- 
eters was formulated. Here, we study this system of equa- 
tions, finding for some important cases analytical expressions 
for local and global order parameters describing synchroniza- 
tion within communities and on the whole network, respec- 
tively. We find that, in the limit of a large number of commu- 
nities, the Ott- Antonsen ansatz introduced in Ref. f23] can be 
used to obtain a low dimensional description of community 
synchrony. Using this description, we characterize the phase 
space of the system where the parameters are the local and 
global coupling. One of our results is that, depending on the 
relative strength of local and global coupling, the transition to 
synchrony in the network can be mediated by local or global 
effects. 

This paper is organized as follows. In Sec. II we describe 
the model. In Sees. Ill and IV we present in detail the local 
and global dimensionality reductions, respectively. In Sec. V 
we discuss the effect of community structure of the network 
on the dynamics and how it promotes hierarchical synchrony. 
In Sec. VI we discuss how our results generalize when certain 
heterogeneities are introduced into the network. In Sec. VII 
we conclude this paper by discussing our results. 



II. MODEL DESCRIPTION 

We are interested in studying coupled oscillators on a net- 
work with strong community structure such that (i) the cou- 
pling strength between oscillators within the same community 
is much larger than the coupling strength between oscillators 
in different communities and (ii) the intrinsic frequency for an 
oscillator is drawn from a distribution specific to the commu- 
nity to which that oscillator belongs. Condition (i) serves as a 
model of situations where all the coupling strengths have sim- 
ilar magnitude, but the density of connections between com- 
munities is less than the density of connections within a com- 
munity. The motivation for condition (ii) is that oscillators in 
different communities could have different frequency distri- 
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butions due to different functional needs (e^-, as in cardiac 
myocytes in different regions of the heart luol ). or as an ap- 
proximation to fluctuations inherent to large but finite com- 
munities. Thus, for a network with C communities labeled 
a = 1,2, ... ,C where community a contains N^- oscillators, 
we assume that the coupling matrix A in Eq. ([T]i can be written 
in block form as Anm = K'^'^ , where a and cr', respectively, 
denote the communities to which oscillators n and m belong. 
Furthermore, we assume that the intrinsic frequencies for os- 
cillators in community a are drawn from a distribution par- 
ticular to that community, denoted by ga{'-^)- We denote the 
fraction of oscillators in community ahy rj„ = N„/N, where 
is the total number of oscillators in the whole network. 

With this notation. E g. (fil l results in the following system, 
considered in Refs. 
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where 9^ denotes the phase of an oscillator in community <t, 
a = 1, . . . C, n ~ 1, . . . , No-, and the intrinsic frequency 
is randomly drawn from the distribution .gcr(w). Next, in order 
to measure synchrony within and between communities we 
define the local and global order parameters 
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Z = i?e'* - 
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respectively, such that measures the degree of local syn- 
chrony in community a and R measures the degree of global 
synchrony over the entire network. We note that the linear 
stability of the incoherent state in this model was studied in 
Ref. il (see also ||2lll ). 



dtfa + de" I fa 
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Following Ott and Antonsen II25I1 . we ex- 
pand fa{0,uj,t) in a Fourier series, f„{0,ui,t) = 

^ (l + Er=i/-,n(^.Oe"' + c.c.), and make the 
ansatz /cr.„(w, t) = a^{uj, t), namely 
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l + ^<Xw,t)e"'' + c.c. , (7) 



which, when introduced in Eq. (|6]l, yields a single ordinary 
differential equation (ODE) 

1 ^ 

ha + lojaa + 2 E (^-'«' - = 0> (8) 

a' = l 

where Za in the continuum limit is given by 
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ga{uj)al{uj,t)duj. (9) 



Finally, by letting the distribution of frequencies ga be 
a Lorentzian with spread Sa and mean fla, i.e. ga{'^) = 
'^cr/{7r[(5^ + (cj — JIct)^]}, we can calculate Za by closing the 
bj contour of integration with the lower-half semicircle of in- 
finite radius in the complex plane and evaluating a* (w, t) at 
the enclosed pole of ga'- 



Za = alifla - iSa,t). 



(10) 



Thus, by evaluating Eq. ^ at u! — fla — iSa, we close the 
dynamics for Za'. 



III. LOCAL DIMENSIONALITY REDUCTION 



In this section, we will study local synchrony by assuming 
there are a large number of oscillators Na in each community. 
Using the definition of Za in Eq. ((Sj, we simplify Eq. Q to 

1 ^ 

(>n - < + ^ E Va'K^^'iza^e-^''^^ - z^.e''''), (5) 

a' = l 

where * denotes complex conjugate. We now move to a con- 
tinuum description by taking the limit N, Na — oo in such 
a way that all rja remain constant. Accordingly, we introduce 
the density function fa{0,uj,t) that represents the density of 
oscillators in community cr with phase 6 and natural frequency 
uj at time t. Since the number of oscillators in each commu- 
nity is conserved, fa satisfies the local continuity equation. 



1 

Za + {Sa - ina)Za + 2 E V'^' i'^' '4 ~ ^'^') = 0' 
a' = l 

(11) 

which defines C complex ODEs, or equivalently 2C real 
ODEs, given by 

Ta = -SaTa + -^"""'Re(^.' e"'"^" ) , (12) 



^a=^a + 4^ y Va'K'"''lmiZa'e-'^'). 

2ra 



(13) 



Equation (fTTT i was formulated originally in Ref. ||25I1 . but its 
consequences for hierarchical synchrony have not been stud- 
ied in detail. Equations ( fT2] i and ( fTsT l describe the dynamics 
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of local synchrony. The synchrony of community a is de- 
scribed by the magnitude of its order parameter r„ and phase 
^a- The phase variable tpa- obeys an equation similar to that 
of the network-coupled Kuramoto model, Eq. ([T]i, but the ef- 
fect of community a' on community a is modulated by the 
degree of synchrony of community a', Va', and its relative 
size ■ In contrast to the Kuramoto model, each community 
has an additional variable which evolves in conjunction 
with the phase variable i/;^- In this sense, the dynamics of the 
community order parameters resembles a network of coupled 
complex Ginzburg-Landau oscillators IztIi . 



In what follows, we will consider the illustrative case in 
which all communities have the same size and spread in natu- 
ral frequencies, i.e. 77^ = C^^ and (5^ = 5- Furthermore, we 
let the coupling strength within each community be the same, 
as well as the coupling strength between oscillators in dif- 
ferent communities. We assume the coupling strength within 
communities is much larger than that between communities, 
namely 



Ck if o- = a' 
K otherwise. 



(14) 



where k and K are of the same order. We clarify that the 
local coupling strength Cfc is chosen so that the local coupling 
within a community is of the same order as the sum of the 
coupling to every other community. More generally, a local 
coupling strength of the form K°'°' ~ k/e with e ^ 1 can 
be analyzed from our results by rescaling k by a factor of Ce. 
In section VI we relax these assumptions and discuss the case 
where community sizes, spread in frequency distributions, and 
coupling strengths vary from community to community. We 
now use the definition of Z in Eq. (01 to rewrite the system in 



State Variables 


Description 




degree of local synchrony of community a 


R 


degree of global synchrony 




phase of local order parameter a 




phase of global order parameter 


Parameters 


Description 


k 


local coupling strength 


K 


global coupling strength 


5 


local frequency spread 


A 


global frequency spread 


0^ 


mean intrinsic frequency of community a 




size of community a 


C 


total number of communities 



TABLE I: Summary of local and global state variables and parame- 
ters of the system. 
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FIG. 1: (Color online) Bifurcation diagram in {K,k) parameter 
space for Eq. (Hi with <5 = A = 1. Re gions A, B, C, and D (described 
in the text) are denoted in red, yellow, green, and blue, respectively, 
with bifurcations (i)-(iv) indicated by solid and dashed curves. 



Eqs. and (O as 
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iij^ = n„ + K 



2r„ 



i?sin(* - -iAcr) 
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We note that although we will let C — > 00 in the next section, 
Eqs. ( fTSl l and ( fT6] l are valid when C is any positive integer 
and can be used to study synchrony on networks with a small 
number of communities. 

Finally, we assume that the mean frequencies $7^ are drawn 
from a distribution G(0), which we assume to be Lorentzian 
with spread A and mean F. However, by entering a rotating 
frame, we can set F = without any loss of generality. For the 
sake of convenience we summarize all local and global state 
variables and parameters of the system in Table U We note 
that choosing a Lorentzian distribution for G{uj) is a natural 
choice if the heterogeneity in the distributions gcr{^) is as- 
sumed to originate from fluctuations arising from the random 
sampling of frequencies from the same Lorentzian distribu- 
tion. In this case, since a sum of Lorentzian random variables 
has a Lorentzian distribution, the distribution of the average 
frequencies in finite communities is Lorentzian. 

Before analyzing Eqs. (fTSl l and ( fTSl l, we illustrate the be- 
havior of the local and global order parameters (5 = A = 
1 over a range of values for K and k. We define r = 

X]ct=i """i^ ^ measure of local synchrony and show the 
behavior of r and R in Fig. [T| While this behavior will be de- 
duced from the analysis that follows, we find it convenient to 
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present the phase space now to provide a framework for our 
subsequent analyses. We note that although the diagram above 
is theoretical, we present plots of R and f following various 
paths in the diagram, and all show excellent agreement with 
the theory. In the parameter space [K, k), we find the follow- 
ing four regions: region A where r,R — (bottom left red), 
region B where r > 0, i? = (top left yellow), and regions C 
and D where r, R > (bottom right green and top right blue, 
respectively). In region A there is neither local nor global 
synchrony, in region B there is local synchrony but no global 
synchrony, and in both regions C and D there is both local and 
global synchrony. We note that although both r,R > in 
both regions C and D, the nature of solutions for are quali- 
tatively different, as will be discussed later Finally, solid and 
dashed curves indicate bifurcations between these regions and 
will be discussed as we proceed with the analysis. In the rest 
of this section, we will study local synchrony, characterized 
by the community order parameters Za- . We will do this by as- 
suming a given value of the global synchrony order parameter 
Z = i?e'*. In the next section, we will study the dynam- 
ics of Z using a dimensionality reduction on the global scale. 
We note here that in the rest of the figures in this paper, since 
we are interested in networks with a large number of commu- 
nities and a large number of oscillators per communities, we 
will compare the results from direct numerical simulation of 
Eq. ^ on networks with large No- and C with the theoretical 
curves obtained from our analysis of the continuum limit. 

First we study local synchrony when i? = 0. In this case, 
from Eqs. (flSl l and ( fTSI l we see that each community decou- 
ples from all others and evolves independently. The phase ipa 
of community a moves with velocity ft^, and the stable fixed 
points of Eq. JTSl l are 
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k~K/C 



if fc - K/C < 25, 
otherwise. 
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so that all are equal. Bifurcation (i), indicated as a solid 
black line in Fig.[Tl is described by this analysis, and occurs at 
k — Kj C = 2(5. To illustrate this bifurcation, we plot in Fig.|2] 
the results of simulating the system as k is varied from zero to 
six with iV^ = C = 400, (5 = A = 1, and fixed = 1 and 
plot the resulting r from simulation (blue circles) against the 
theoretical prediction of Eq. ( [TtI i (dashed red). The interpre- 
tation of this result is that the oscillators in each community 
synchronize as in the all-to-all Kuramoto model, but with an 
effective coupling strength k — K/C, which shows that the 
weak coupling to other independently evolving communities 
slightly inhibits synchrony. 

The analysis above assumes i? = 0. Now, we will ana- 
lyze local synchrony when i? > 0. In this case some of the 
communities become synchronized with each other Given a 
value of Z (which can be obtained using another dimension- 
ality reduction, as we will show in the next section), com- 
munity a synchronizes with the mean field [i.e., a solution 
■0(7 = 0, Tcr = for Eqs. ( fTsl l and (fTSl i exists] if 
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FIG. 2: (Color online) Average degree of local synchrony r versus k 
from simulation (blue circles) with C = A^o- = 400, (5 = A = 1 and 
K = 1 compared to theoretical prediction in Eq. l llVb (dashed red). 



in which case 



ip„ — ~ arcsin 



KR{rl 
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(19) 



and otherwise the community drifts indefinitely. The degree 
of local synchrony for locked communities can be found 
by setting ?'> in Eq. ( fTSl l to zero and using Eq. ( fT9] l, which 
gives the implicit equation 
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Eq. (|20] i determines the steady-state value of for locked 
communities and yields two possible kinds of solutions for 
r^: either Eq. ( l2Qt has a real solution for every 11^, or it 
has a real solution for only some It can be shown that 
when k — K/C < 26, Eq. ( l20t has a real solution for all 
ila, and thus each community becomes phase-locked and each 

reaches a fixed point as t — > oo. On the other hand, if 
k — K/ C > 26, there is a real solution for only some fl^r 
with magnitude less than a critical locking frequency, which 
we denote as il. In this case communities with {fla-l < 
phase-lock and is given by the solution of Eq. ( |20l ). while 
other communities continue drifting indefinitely. The phase 
angle i/'o- of a drifting community a increases or decreases 
mono tonic ally and therefore its order parameter might be 
time dependent, according to Eq. ( fTSl ). However, assuming a 
stationary global order parameter with constant R and 5* (as 
will be discussed in the next section), the solution of the two- 
dimensional autonomous system in Eqs. ( fTSl ) and ( fTSI ) must 
approach a limit cycle (this can be shown, for example, us- 
ing the Poincare-Bendixson theorem ll28ll ). To estimate the 
time averaged value of in this limit cycle, we neglect the 
effect of the cosine term in Eq. ( fTSb over one period and find 
that the time averaged order parameter for drifting commu- 
nities is approximated by {r^) ~ J\ — k~^K/c • ^^^^ value 
agrees with the solution of Eq. ( [20l i when 0,^ is the locking 
frequency in Eq. (fTSl l. Therefore, the community locking fre- 
quency can be determined by inserting the expression for 
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FIG. 3: (Color online) Time- averaged r vs Q from simulation (blue 
circles) with C = = 400 and 5 — A — 1 compared to theoret- 
ical prediction (dashed red) for {K, k) — (1, 8) (a), (8, 1) (b), and 
(4, 8) (c). The vertical arrows indicate the theoretical value for the 
locking frequency obtained from Eq. l l21b . 



above into Eq. ( fTSl l. obtaining that communities lock when 
their frequency il„ satisfies 
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The locking frequency only is defined for k — K/ C > 25, 
which defines a new bifurcation. When k—K/C > 2S (region 
D), the locking frequency fl is finite and only some communi- 
ties phase-lock. As fc — K/C approaches 26 from above, the 
locking frequency diverges. For fc — K/C < 2d (region C), 
all communities phase lock. The boundary between these two 
regions for larger K is denoted as bifurcation (ii), and is indi- 
cated as a solid black line in Fig. 1. A heuristic interpretation 
of this transition is that when fc is increased through bifurca- 
tion (ii), communities with large \ desynchronize because 
the local coupling strength k causes them to prefer an angular 
velocity ^ much closer to their own mean frequency fl^ than 
the mean frequency of the entire network. 

To test Eqs. ( fTTl ), ( |20l i. and ( |2TI ). we simulate the system 
with = C = 400 and (5 = A = 1 with {K, k) = (1, 8), 
(8, 1), and (4, 8) (parameters from regions B, C, and D, re- 
spectively) and plot time-averaged as a function of Vt^r in 
Figs. [3|a), (b), and (c), respectively. Results from direct sim- 
ulation are plotted in blue circles and compared to theoretical 



predictions, which are plotted as dashed red curves. Fig.[3ta) 
corresponds to region B, where R ~ Q and r„ is given by 
Eq. ( [TtI i and is therefore independent of a. Fig. [3|b) cor- 
responds to region C, where all communities lock and their 
order parameter is a solution of Eq. ( l20b . Fig.[3jc) corre- 
sponds to region D, where some communities lock and their 
order parameter Tr, is a solution of Eq. ( [2Qb , and other com- 
munities drift and their order parameter is independent of 
(T and given by (7^). The vertical arrows indicate the theoret- 
ical value for the locking frequency obtained from Eq. ( |2TI ). 
Theoretical results match very well with the numerical simu- 
lations. 



IV. GLOBAL DIMENSIONALITY REDUCTION 



In the previous section, we studied local synchrony by as- 
suming a steady-state value for the global synchrony order pa- 
rameter Z = i?e**. We now discuss how the global order 
parameter can be found by making a second dimensionality 
reduction on a global scale. As we previously let N^r tend 
to infinity in order to enter a continuum description within 
each community, we now consider the limit C — > 00 and in- 
troduce the density function fi, r, t) that describes the 
density of communities with average phase ip, mean natural 
frequency il, and degree of local synchrony r at time t. In 
analogy with individual oscillators, the number of communi- 
ties is conserved and F must satisfy the continuity equation 
dtF + d^,{F^) + dr{Fr) = 0. However, we find that the 
degrees of local synchrony r quickly reach a stationary distri- 
bution, so we seek solutions where dr{Ff) = 0. In analogy 
to the classical Kuramoto model, we find that approaches 
a fixed point if community a phase-locks, or otherwise forms 
a stationary distribution with other drifting r's. With Eq. ( fT6l ) 
the continuity equation becomes 



dtF + d^\F 



n + K 



2r J 



-] lm{Ze-'"f') 



0, (22) 



where r = r{il, R) is the steady-state value of r given by 
Eq. {TT} or implicitly by Eq. ( |20l ). 

Like Eq. (|6]l, E q. (l22b is of the form studied by Ott and 
Antonsen in Refs. 125112911 . and can be solved with a similar 
ansatz. Thus, we make the ansatz 



F{ib,n,r,t) 



2tt 



l + ^A"(l],r,t)e™''' + C.C. 



(23) 



Inserting Eq. ( |23] | into Eq. ( |22] |. we find that 

K fr'^ + V 



A + iVlA 



4 V r 



{A^Z-Z*)=Q. (24) 
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FIG. 4: (Color online) Degree of global synchrony R (main) and 
average local synchrony r (inset) versus K from simulation (blue 
circles) with TV^ = C — 400, 5 — A = 1, and k — 4 compared 
to theoretical prediction from Eqs. ( I30t . Bib . ( I32t and B3b (dashed 
red). 



We calculate Z as: 
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Gin)A*{n,r,t)rdn. 
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Since r(r2, i?) is defined implicitly by Eq. ( l20l i for locked 
communities and by Eq. ( fTTI ) for any drifting communities, 
it is potentially piecewise-defined and not smooth. However, 
to a very good approximation we can do this integral using 
residues by considering the solution r(r2, R) of Eq. ( |20| ) for 
< which is real and positive for Im($7) 0~ as a func- 
tion of complex fl. The function r is analytic when Im(il) < 
0, and its real part converges to r{il, R) as lm(f7) 0^ with 
|0| < O, while its imaginary part converges to an odd func- 
tion. As lm{fl) 0~ for \Q\ > h, the real part of F dif- 
fers from Eq. (fTTb by a bounded amount. If G{il) decays so 
quickly that the error in approximating r by r for > il can 
be neglected when computing the integral, we can approxi- 
mate the integral above by the integral which has r instead 
of r (due to the symmetry of G, the imaginary part of r does 
not contribute to the integral). The integral with r on the real 
line can be done by deforming the contour of integration to 
the line connecting zi = — B — ie to Z2 = B — ie, where B, 
e > 0, and closing the contour with the semi-circle in the neg- 
ative complex plane connecting Z2 to zi. Using the residue 
theorem, and taking i? — > oo and e 0, we obtain 



Z K?A*{-iA,?,t), 



(26) 



where we have defined r = ?■(— iA, R). For the Lorenzian 
distribution with A = 1, we expect this approximation to be 
excellent when il > 4, but the agreement between the direct 
numerical simulation of Eqs. (|2]i and the theoretical predic- 
tions is very good even for situations in which fl is smaller 
[e.g., Fig.|5](a) close to the transition for R]. We note that if 
{k, K) is in region C [see Fig. [T| using r to evaluate the in- 
tegral in Eq. jZSl ) is exact since all communities lock and are 
described by Eq. ( l20l i. 



Evaluating Eq. ( |24] | at 
complex dynamics for Z: 



-iA and r = r closes the 
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The evolution of R and 5* are given by 



Ai?+^i?(f^ + l) 



(27) 

(28) 
(29) 



We note that these equations are valid provided that (a) F is in 
the manifold of Poisson kernels [i.e. is of the form in Eq. (l23l ll 
and (b) the distribution of degrees of local synchrony r re- 
mains stationary as the system evolves. Regarding assumption 
(a), Ref. ll29ll shows that in the Kuramoto model all solutions 
approach this manifold as t — > oo. The stable fixed points of 
Eq. (|28l) are 
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otherwise. 
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To elimina te r we assume nonzero R (and thus r > 
^4 A/ A' - 1), and insert Eq. ^ into Eq. ^ with fla = 
— lA. We choose the real, positive solution given by 



A-S+ ^{k + K-5Y - 2{k + K + (5) A + A2 

'^^ V VTk ' 

(31) 

which we insert back into Eq. (|30] l to obtain R. We note that 
other solutions for r are purely imaginary or negative. From 
the top line of Eq. ( l30l l, the imaginary solutions for r result 
in a critical value for K larger than 4A, while real solutions 
result in a critical value smaller than 4A, and thus we choose 
the positive real solution (the negative solution results in i? < 
0). Finally, to calculate the bifurcation curve for the onset of 

global synchrony, we let r 4A jK — 1 which yields 

the curve fc = k^is ^ t- ''"^'^ curve is indicated as a dashed 
black curve in Fig. [T]and gives bifurcation (iii) from region A 
to C and bifurcation (iv) from region B to D. 

We now seek to compute the mean degree of local syn- 
chrony r. In the large C limit we consider here, r is given 
by an integral equation. If (A', fc) is in region C, i.e. fc < 25, 
then since each community becomes phase-locked, we simply 
have 



G(Sl)r{9.,R)d9.. 



(32) 



However, if (A', fc) is in region D, i.e. k > 2S, then because 
some communities phase lock and some do not, we have that 



r = ; ^ G{n)r{n,R)dn 
'\n\<n{R) 



G{n){r)dQ, 
(33) 
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FIG. 5: (Color online) Degrees of global synchrony R (blue circles) 
and average local synchrony r (red triangle) along paths (a) k — 
3K/2 and (b) k = K/2 from simulation with = C = 1000 and 
5 = A = 1. 



where il is the locking frequency given by Eq. ( |2TI ). 

To illustrate these results, we simulate the system with 
Na ^ C ^ 400, 5 = A = 1, fc = 4, and let K vary be- 
tween zero and six. In Fig. 2] we plot R (main) and f (inset) 
from simulation in blue circles and the theoretical predictions 
from Eqs. ( l30l l. (ISTT i. ( |32] | and (l33T l in dashed red. Theoretical 
predictions agree well with simulations. 



ical predictions plotted in black dashed and dot-dashed, re- 
spectively. Note that for these parameters = 1/ V^, so 
mi > Trie > and accordingly we see a separation of local 
and global onset in Fig.|5ja), but not in Fig.|5jb). 

We interpret these results as follows. Along paths where 
k > rricK local coupling effects dominate global coupling 
effects. In this case the community structure is strong enough 
to yield a hierarchical ordering of synchrony, i.e. a separation 
in the onset of local and global synchrony. However, when 
k < rricK global coupling effects dominate local coupling 
effects. In this case the community structure is weak enough 
to yield a simultaneous onset of local and global synchrony. 



VI. HETEROGENEITIES 

We now discuss how the results above generalize when 
some of the assumptions previously used are relaxed. We al- 
low for heterogeneities in both the sizes of communities and 
spread in frequency distributions gcr(w), i.e. we allow 77^ and 
5a to vary from community to community. We also allow the 
local and global coupling strengths to vary, letting K'^'^ = k'^ 
for a = a' and K'^ for a ^ a' . 

Beginning with local synchrony, we carry out a dimension- 
ality reduction on the local scale and obtain the following 
ODEs: 

/ K^X 1-7-2 

ra = ~5ara + Crj^ Ik" — j — ^ 



^a^na + 



2r, 



(34) 

i?sin(*-V^<,). (35) 



V. HIERARCHICAL SYNCHRONY 

With a complete understanding of both local and global 
synchrony in the system studied above, we now discuss hi- 
erarchical synchrony. We consider moving slowly (compared 
with A^^) along some path in {K,k) parameter space, re- 
stricting paths to lines starting at (0, 0) for simplicity. From 
our analysis we find that bifurcations intersect at {K, k) = 
(A - <5 + VA2 + + 6A(5, 25). Thus, for Hnes k = mK, if 



m > rrir 



25 



- — ; — , ^ the onset of local synchrony 
occurs before the onset of global synchrony. On the other 
hand, if tti < nic the onset of local and global synchrony oc- 
cur simultaneously. Choosing mi = 3/2 and ?Ti2 = 1/2, 
in Figs. |5] (a) and (b) we plot the steady-state values of R 
and r resulting from moving along the lines k ~ miK and 
k = m2K, respectively, for ^ C = 1000 and (5 = A = 1. 
We note that = C = 1000 is used in these simulations 
rather than 400 as in the previous simulations because we 
find that finite-size effects are more prevalent near bifurcation 
(iii). This is most likely due to the fact that at this bifurca- 
tion the onset of local and global synchrony occurs simulta- 
neously. The values of R and r from simulation are plotted 
in blue circles and red triangles, respectively, with theoret- 



2S^ 

Cr^^ik'-K" /C) 



Thus, when i? = 0, we have that 

if Cr]a{k'^ - K^/C) < 25a 

otherwise. 

(36) 

The onset of local synchrony in community a occurs at 
k'^ — K°' /C = 25a-/ {Crja), i-e- in general synchrony occurs 
at different values for different communities. When i? > 0, 
community a becomes phase-locked if 



+ 1 

10 I < K^R^ 



in which case r„ satisfies 



2ra 



1 -2 



(37) 



daVa = Cl^a [ k — I Ta 



K"R- 



K'^'^R^rl + iy 



(38) 



otherwise community <t will drift. Note that for a given value 
of R, the behavior of Va depends not only on ft a, but also 
Cr]a, 5a, k"', and A''^, so in general there is no single locking 
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frequency that separates locked and drifting communities at 

To study global synchrony, we again perform a dimension- 
ality reduction on the global scale. Since rja-, 5^,, k"', and K"' 
vary from community to community, after sending C — > oo 
we introduce the density function F^i/j, fl, r, rj, S, k, K, t) that 
represents the fraction of communities with phase ip, mean 
natural frequency il, degree of local synchrony r, size r], fre- 
quency distribution spread 5, and local and global coupling 
strengths k and K at time t. Noting that Va depends on rja, 
5cr, k"', and K'^ and again looking for solutions with station- 
ary Tcr, F satisfies the continuity equation 



dtF + d^iF 



n + K 



1 



2r 



lm(Ze 



= 0, (39) 



where now r depends on 77, 6, k, and K in addition to fl and R. 
We assume that for each community the mean frequency O^, 
size frequency distribution spread S^, and local and global 
coupling strengths k°' and K"^ are all chosen independently 
and make the ansatz 



Fi4,,n,r,7j,d,k,K,t) 



G{n)H{r])D{6)J{k)L{K) 
2^ 



1 + J2 ^"(^' ^' ^' + c-c. , (40) 



which yields the ODE 



A + iflA- 



K / + 1 



{A^Z - Z*) = 0. (41) 



Finally in the continuum Umit Z can be calculated by the 
integral 

/•oo /"OO /"OO />! /"OO /'27r 

Z= / / / / / F{i^,n,r,Tj,6,k,K,t) 

Jo Jo Jo Jo J -00 Jo 

X re'''' dijdfldrjdSdkdK 

00 poo nOO nl 

/ / / HirmS)J{k)L{K) 
Jo Jo Jo 

X ?A* (-ir, f , T], S, fc, K, t)dr]dddkdK. (42) 

Equations. i4li and ( l42l i govern the global synchrony of 
the system and must be solved self-consistently with the local 
dynamics, governed by Eqs. ( |34] | and (l35l l. For arbitrary dis- 
tribution functions H{rj), D{S), J{k), and L{K) the integral 
in Eq. (l42l i might need to be evaluated numerically, but for 



certain choices, e.g. exponentials or linear combinations of 
Dirac delta functions, further analytical results are attainable 
but not presented here. 



VII. DISCUSSION 

We have described and solved fully the steady-state dynam- 
ics of coupled phase oscillators on a modular network with a 
large number of oscillators in each community and a large 
number of communities. In particular, we have studied lo- 
cal and global synchrony, i.e. synchrony within and between 
communities, respectively. First we assumed a large num- 
ber of oscillators in each community and used a local dimen- 
sionality reduction to study local synchrony. Next, when the 
number of communities is large, we showed that a global di- 
mensionality reduction can be done to study global synchrony. 
Our analytical results shed light on the phenomenon of hierar- 
chical synchrony, characterized by synchronization on a local 
scale before it occurs on a global scale, which occurs when 
the community structure of the network is strong enough. The 
system analyzed in this paper modeled synchrony on a net- 
work with two community levels, but synchrony on networks 
with more levels, e.g. communities with subcommunities, can 
be modeled in a similar way and analogous analytical results 
can be obtained. 

Although we have assumed strong uniform coupling within 
communities and weak uniform coupling between communi- 
ties, we conjecture, based on preliminary numerical experi- 
ments, that the system studied in this paper is in some cases a 
good quantitative model for networks where links between os- 
cillators in the same community are dense and links between 
oscillators in different communities are sparse. 

An interesting result is that the system of planar oscillators 
representing community interactions [Eqs. (fTsT i and ( fT6] l1 ad- 
mits an approximate low dimensional description. The anal- 
ysis of community synchrony in Sec. IV is, to the best of our 
knowledge, the first low-dimensional description of oscillator 
systems in which each oscillator has a phase and an associ- 
ated oscillation amplitude. Other systems of coupled planar 
oscillators could be analyzed in the same way. 
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